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Abstract: Magnetically-doped graphene systems are potential candidates for application 

>- , in future spintronic devices. A key step is to understand the pairwise interactions between 

^ . magnetic impurities embedded in graphene that are mediated by the graphene conduction 

^ ■ electrons. A large number of studies have been undertaken to investigate the indirect 

exchange, or RKKY (Ruderman-Kittel-Kasuya-Yosida), interactions in graphene. Many of 

O ■ these studies report a decay rate faster than expected for a two-dimensional material and 

the absence of the usual distance dependent oscillations. In this review we summarize the 

techniques used to calculate the interaction and present the key results obtained to date. 

X : The effects of more detailed parameterisations of the magnetic impurities and graphene host 

, 

o3 . are considered, as are results obtained from ab initio calculations. Since the fast decay of 

the interaction presents an obstacle to spintronic applications, we focus in particular on the 
possibility of augmenting the interaction range by a number of methods including doping, 
spin precession and the application of strain. 

Keywords: graphene; exchange interactions; RKKY; spintronics; magnetic impurities; 
magnetic interactions 
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1. Introduction 

Graphene has been in the scientific limelight for the past few years due to a range of interesting 
properties that it has been shown, or predicted, to display. Superlative mechanical properties and a unique 
set of potentially tunable electronic and optical properties suggest a large number of possible applications 
for graphene and related materials [1-4]. Apart from potential roles in nanoelectronics, optoelectronics 
and as reinforcements in composites, graphene-based materials have also been mooted for application 
in spintronics [5]. This field is projected to play a very important role in future technologies as it 
may allow information storage, processing and communication at faster speeds, and with lower energy 
consumption, than is possible with conventional electronics [6-9]. The principal idea of spintronics is 
to exploit not only the charge of an electron, but also the spin degree of freedom associated with its 
intrinsic angular momentum. Despite carbon not being magnetic, graphene-based spintronics may be 
achievable when we consider that many of its derivative materials and nanostructures display various 
scenarios of magnetism [5]. Graphene has also many properties that suggest its possible use as a carrier 
of spin information, including weak spin-orbit and hyperfine couplings, which are the main sources of 
relaxation and decoherence of electron spin [10-16]. 

Many of the proposed graphene-based spintronic devices are underpinned by the spin-polarised edge 
state that is predicted to occur when a graphene sheet is cut to have the so-called zigzag edge geometry. 
Particular focus has been paid to zigzag-edged graphene nanoribbons, narrow strips of graphene with 
parallel zigzag edges that are predicted to display opposite spin orientations. Such a system may allow 
a possibility of tuning its spin-transport properties [17-22]. For example, the prospect of triggering a 
half-metallic state using extemal electric fields in zigzag-edged nanoribbons has been suggested [19]. 
The realisation of such a device would allow efficient electronic control of spin transport, a very useful 
property in spintronics and something that is difficult to achieve in other materials. Despite theoretical 
advances in the study of nanoribbons, and some recent experimental results hinting at signatures of this 
edge magnetism [23], the difficulty in patteming the edge geometries required for these effects to be 
observed may prevent a wider scale exploitation. Furthermore, the spin-polarised edge state for zigzag 
edges is predicted to be highly dependent on the edge geometry and not particularly robust under the 
introduction of edge disorder [24]. These factors present major obstacles in the path of utilising the 
intrinsic magnetic edge states of graphene in experimentally realisable devices. 

Another possibility that has been proposed for graphene-based spintronics is the exploitation of 
defect-driven magnetic moments that arise in graphene [25-27]. Magnetic moments have been predicted 
to form around vacancies and other defects in the graphene lattice and the possibility of engineering a 
ferromagnetic state in graphene from such moments has been suggested. However, such a claim would 
seem to be restricted by the implications of the Lieb theorem [28], which states that any such magnetic 
moments arise from a disparity between the two sublattices of graphene. Large-scale, randomised 
disorder would tend to minimise such a disparity and prevent the formation of a ferromagnetic state. 
However, recent experimental evidence suggests the possibly of engineering such a state through partial 
hydrogenation [29]. The existence of such a state may then be accessible through magnetoresistance 
measurements [? ]. 
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A third possibility for incorporating graphene in spintronic devices, and one that we shall focus on in 
this review, is through the doping of graphene with magnetic impurity atoms, as shown in Figure 1. 
This approach allows graphene to act as host to an ensemble of transition metal atoms, or indeed 
more generic magnetic species, and mediate the interactions between them, potentially allowing for 
long ranged ordering or the transfer of spin information between impurities. In order to predict the 
magnetic properties of such systems, it is essential to understand both how a single impurity connects to 
the graphene lattice and also the nature of the interactions between multiple impurities. To address the 
former, a large number of studies have investigated the hybridisation and resultant electronic, magnetic 
and magnetotransport properties of graphene systems with different magnetic species embedded in them 
[31-33,35? ]. However, the interaction between the localised moments formed at the impurity sites 
is predicted to be largely independent of the magnetic impurity species chosen and instead depend 
strongly on the electronic structure of the host material. This is because the indirect exchange coupling 
(lEC) [36-39], often referred to as the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction [40-44], 
between magnetic impurities in a graphene system is mediated by the conduction electrons of the 
graphene host. Throughout this work we will use the terms "lEC" and "RKKY" interchangeably to 
refer to a general conduction-electron mediated interaction between magnetic objects, unless specifically 
referring to the RKKY approximation to the interaction. A long-ranged interaction of this type allows 
impurity moments on graphene to feel each other's presence and respond to magnetic perturbations or 
excitations at other impurity sites. The usual manifestation of this interaction is through the lowering 
of the total energy of the system when the moments adopt certain favourable alignments. In this way 
the interaction can dictate any magnetic ordering that arises between the moments. These types of 
interactions have been investigated previously in multilayer devices, leading to the discovery of the Giant 
Magnetoresistance effect [45,46], and more recently in materials closely related to graphene, namely 
carbon nanotubes [47-49] . 

To exploit the interaction between impurities, it is important to understand its behaviour as the 
separation between them is varied. Previous studies suggest that this type of interaction should 
oscillate and decay as the impurity separation is increased [36-44]. Furthermore, the oscillation 
period is determined by the Fermi surface of the host material, whereas the decay rate depends on 
the dimensionality of the host and magnetic objects. For impurity atoms in a one-dimensional or 
two-dimensional host medium, the interaction should be expected to decay as or respectively, 
where D is the distance between the magnetic impurities. The decay rate in one-dimensional systems 
is borne out by investigations in metallic carbon nanotubes where the predicted long range decay is 
found theoretically for substitutional magnetic impurities [47]. However, unique features arising from 
the peculiar electronic structure of the underlying graphene lattice were also reported. The sign of the 
interaction, determining the preferred alignment of the moments on two impurities, was found not to 
oscillate as a function of distance as seen in other materials, but rather to depend only on the connections 
of the impurities to the two triangular sublattices composing graphene. These sublattices are represented 
by filled and hollow circles in Figure 1 . Furthermore, for certain impurity configurations, the interaction 
was found to decay much faster than expected [48-50]. 

Similar features to those reported for carbon nanotubes have now been predicted in graphene. Due 
to the interest in carbon-based spintronics, a large number of studies in recent years have investigated 
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indirect exchange interactions between localised moments in graphene. Sublattice-dependent signs and 
faster than expected decay rates have been seen when a variety of methods are employed, ranging from 
analytic or numerical treatments of the interaction using tight-binding or Dirac Hamiltonians to more 
intensive total energy calculations within a density functional theory (DFT) framework. In this review 
the details of the various approaches will be discussed in Section 2 and the consensus behaviour for 
simple impurities that has emerged from the majority of the studies analysed in Section 3. The results 
for impurities beyond the simplest substitutional case and for considerations often neglected by the basic 
RKKY approximation are discussed in Section 4. Finally, in Section 5 we consider how the interaction 
between magnetic impurities can be manipulated in order to extend the interaction range or otherwise 
alter its behaviour. 

Figure 1. Schematic illustration of two substitutional magnetic impurities (red circles) 
embedded in a graphene lattice. Also shown are the high symmetry armchair (A) and zigzag 
(Z) directions, and the two-atom unit cell and lattice vectors ai and a2. Atoms from the two 
triangular sublattices are represented by filled and hollow symbols. 




I \ / 



►A 

2. Methods 

In this section we will outline several of the general methodologies employed to calculate the indirect 
exchange coupling in graphene. We begin by considering an energy difference calculation between the 
parallel and antiparallel orientations of the localised moments on two magnetic impurities located at 
sites A and B in the graphene lattice. These configurations are represented schematically in Figure 2. 
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Employing the Lloyd formula [51] to calculate the total energy difference between these two spin 
configurations we arrive at the Quantum Well method result proposed by Edwards et al. [38,39]. We 
shall briefly derive this result in a completely generic form in terms of Green functions without stating the 
form of the Hamiltonian chosen to describe the electronic structure of the graphene host or the magnetic 
impurities. From this result we identify the common approximations often made in the literature and 
demonstrate how the commonly used RKKY expression can be reached. 

From these general results we turn our attention to the specific case of graphene. Two common 
Hamiltonians used to describe the electronic structure, namely the tight-binding and linearised Dirac 
formalisms, are introduced and compared. The various approaches, numerical and analytical, to 
calculating the coupling within these formalisms are discussed. Finally, the merits and disadvantages 
of these methods, and those of more complete ab initio treatments, are given. 

Figure 2. Schematic representation of the spin-split potentials at the magnetic sites for the 
parallel (a) and antiparallel (b) alignments of the magnetic moments. The red (blue) curves 
show the potential experienced by up-spin (down-spin) electrons as a function of position. 
For simplicity the host medium is represented by a one-dimensional linear chain. 




b) 
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2.1. Calculating Exchange Interactions 

The lEC between a pair of magnetic objects is usually defined as the energy difference between the 
parallel and antiparallel alignments of the localised magnetic moments on the objects. We distinguish 
the indirect interaction by assuming that it is mediated entirely by the conduction electrons of the host 
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material, in contrast to the direct exchange which depends on the overlap of the impurity electron 
orbitals and which decays very abruptly. To begin, we consider an initial configuration of two magnetic 
impurities at sites A and B, separated by a distance D, in the host medium. We introduce a Green 
function, g, that describes the electronic structure of this system. For the moment we will keep the 
form of this Green function, and its associated Hamiltonian, as general as possible. We assume that it 
describes in sufficient detail electron propagation in the host material and also spin-dependent potentials 
localised at the magnetic impurities, which reflect the fact that electrons in the system are subject to a 
different potential in the magnetic regions than elsewhere, as shown schematically in Figure 2. The 
resultant potential profiles seen by "up" and "down" spin electrons are the result of solving the Hubbard 
Hamiltonian and are illustrated by red and blue lines respectively. The allowed energies of the system 
are quantised and furthermore, changing the separation between the magnetic objects shifts the allowed 
energies relative to the Fermi energy of the system giving rise to an oscillatory coupling [52]. The 
coupling can be calculated by summing over all the energy levels below the Fermi energy and taking 
the difference between the cases with spin potentials corresponding to moments aligned parallel or 
anti-parallel. In this manner, the calculation of the lEC reduces to the calculation of an energy difference 
between two distinct configurations. Such a calculation can be simplified using the Lloyd formula [51], 
which allows the calculation of the energy difference directly without the need to calculate the total 
energy of either configuration. 

The Lloyd formula expression to calculate the total change in the energy of a system due to a 
perturbation is given by 

AE{EF) = ^lra J dE f{E) In (det{i - g{E)V)^ (1) 

where g as above is the Green function describing the unperturbed system, V is the applied perturbation 
potential and f{E) is the Fermi function. The unperturbed system consists of two moments embedded 
at sites A and B, and aligned parallel along the z-direction so that the angle between them, ^ = 0. This 
setup is shown schematically in the panel a) of Figure 2. For simplicity we consider one magnetic orbital 
on the impurity sites, but the calculation can be easily generalised to a more realistic magnetic impurity. 
We associate a bandcentre (6) and exchange splitting (Vex) with this orbital as illustrated in Figure 2. 

We now introduce a spin perturbation which rotates the magnetic moment at B by an angle 9 with 
respect to that at A. This perturbation is given by 



V{9) = -Vex [(cos^- l)<T^ + sin^a^ 



cos 6 — 1 sin 6 
sin 6 1 — cos 6 



(2) 



where &z and are the relevant Pauli matrices. Since the magnetic moments break the spin degeneracy 
of the electrons, the Green functions must be written in terms of their up- and down-spin components. In 
the initial coUinear configuration there is no mixing between the spin bands and (^bb is diagonal in spin 
space. It can be shown that 



det 



/ - gsBiE) vie)] = 1 - 2Vl gU 9^3 (cos 6 - 1) (3) 
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where g\g ((7^^) are the off-diagonal elements of the up (down) spin Green function connecting sites 
A and B. Taking = tx, corresponding to an antiparallel alignment of the moments we find, for zero 
temperature, that 



1 r / 

Jba = -AE{e = n) = —lm In 1 + 4 gl^{E) g\^{E) 



(4) 



where we have chosen the sign convention where a negative value of the coupling, J, corresponds 
to a preferential parallel alignment of the moments and a positive value to a preferential antiparallel 
alignment. When performing the integral in Equation (4) numerically we can rewrite the integral over 
the imaginary axis where the integrand tends to be smoother and easier to integrate. With an imaginary 
axis integration, the expression for the coupling becomes 



1 1"°° 

Jba = - In 1 + 4 gl^iEp + ly) gieiEp + ly) 

Jr, 



(5) 



We note that in these expressions the distance dependence is entirely contained within the product of 
off-diagonal Greens function elements and the behaviour of these quantities will dictate the behaviour of 
the interaction as the separation between the moments is varied. 

The expression for the coupling given in Equation (4), with the log term in the integrand, is exact. That 
is, it gives the exact energy difference between the EM and AEM configurations of the system described 
by the Green functions used. In the next section we will discuss how the relevant Green functions will be 
calculated, but first we focus on a common perturbative approach to calculating the coupling. To proceed 
analytically it is useful to note that Equation (4) can be written as a perturbation expansion in powers of 
the exchange splitting Vex and when expressed to leading order in V^x gives 

Jba = -^ j dEf{E)\m[gBA{E)gAB{E)] (6) 

where gsAiE) is the spin-independent Green Eunction describing electron propagation in the pristine 
host material [53]. 

This expression for the coupling is equivalent to the Ruderman-Kittel-Kasuya-Yosida (RKKY) 
approach, initially developed to describe the coupling mechanism of nuclear magnetic moments [40], 
then expanded to describe a wider range of indirect coupling phenomena [41,42] and generalised to 
provide a model capable of reproducing experimental observations [43,44]. This treats the coupling as 
a consequence of the spin polarisation of the conduction electrons of the host by the magnetic objects. 
Under this approach, the coupling is written as an effective direct coupling, Jba sb ■ s^, between the 
two moments, where the coupling strength is given by 

Jba = (^) xIa (7) 

where A is an adjustable parameter representing the magnitude of the coupling between localised spins 
and conduction electrons and x is the static magnetic susceptibility which relates the response of the host 
magnetisation to a static magnetic field. Writing this quantity in terms of Green functions, 

Xba = -1 J dE f{E)lm[gBA{E)gAB{E)] (8) 



Crystals 2013, 3 



8 



and the equivalence of Equations (6) and (7) becomes clear. 

Thus the RKKY approach can be considered a second-order perturbational approximation to the lEC. 
The approaches are in general equivalent, except when Vex >> 1 or when the host material possesses 
certain symmetries [53-56]. The RKKY approximation generally provides a good description at large 
separations when the coupling is quite small. An important contrast between the Lloyd formula approach 
and the RKKY approach is the choice of Green function used. The former approach avails of the 
spin dependent Green functions in the ferromagnetic configuration, g\j^ and g\g whereas the latter 
approach generally uses their pristine, spin-independent counterparts. We note that use of the pristine 
GFs essentially sets the value of the bandcentre, 5, to zero. An intermediate approach can be taken by 
using the spin-dependent Green functions with an RKKY-like expansion of the logarithm in the Lloyd 
expression. The general form expected for the coupling is 

■HD) ^ (9, 

where kp is the Fermi wavevector in the conducting host and the exponent a determining the rate of 
decay with separation depends on the dimensionality of the system. 

Most of the analytic investigations into the lEC in graphene make use of the RKKY approximation. 
Although we have written it here in terms of energy dependent Green functions, it is often formulated 
differently. A common alternative formulation is in terms of time dependent Green functions, the Fourier 
transforms of the GFs used here, although conflicting results have been reported depending on whether 
the real or imaginary time representations are used [57-59]. Numerical approaches to indirect exchange 
interaction generally proceed in one of two ways. The first uses numerical methods to calculate either 
perturbative or exact expressions for the coupling like those given above. Another method involves direct 
calculation of the energies of the individual FM and AFM configurations. This is a common approach 
in ab initio calculations [60-64] although it has been employed also for tight-binding models [65,66]. 
The computational cost of this type of calculation tends to limit these investigations to small or medium 
values of separation and it can be difficult to perform a meaningful comparison with results for the 
asymptotic behaviour calculated in other ways. 

2.2. Calculations in Graphene 

Now that we have introduced the general methodology to calculate indirect exchange couplings, we 
shift our focus to examine the specific case of interactions in graphene. The electronic structure of 
graphene plays a key role in determining the interaction properties and so we first outline the different 
methods used to calculate or approximate the band structure. We then discuss how the Green functions 
required for some coupling expressions can be found and compare the different approaches used to 
calculate the coupling. 

The electronic structure of graphene near the Fermi energy is principally determined by tt bonds 
formed from the out of plane 2pz orbitals of carbon atoms situated on the vertices of a hexagonal, 
or honeycomb, lattice. The remaining carbon orbitals hybridize in plane to form very stable a bonds 
which give graphene its extraordinary strength. Thus a single-orbital nearest-neighbour tight-binding 
Hamiltonian is usually more than adequate when describing the electronic properties of graphene. To 
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write this Hamiltonian we make use of the two-atom unit cell shown in Figure 1 . Thus every site in the 
lattice is denoted by a vector r giving the unit cell location and an index n = 0,1 differentiating between 
the two sites within the unit cell, illustrated by filled and hollow symbols in the figure. The tight-binding 
Hamiltonian is written 

H= ^ |m)t::;'(rV| (10) 

<r,n,r',n'> 

where |rn) represents the orbital at the site given by {r,n} and t^'^, is the electronic hopping term 
between two such orbitals. The sum is restricted to orbitals at neighbouring sites and the non-zero 
hopping terms all take the value t = — 2.7eV. One consequence of only allowing nearest neighbour 
hoppings is that graphene is a bipartite lattice. The sets of filled and hollow symbols in Figure 1 can 
be regarded as intersecting triangular sublattices that together form the honeycomb graphene lattice. 
Introducing hopping terms beyond the first nearest neighbour break the bipartiteness of graphene. 
However, since these terms are much smaller than the nearest neighbour hoppings, graphene behaves 
effectively as a bipartite lattice. It is this feature of the graphene lattice, combined with the half filling of 
its TT orbitals, that determine many of the features of the RKKY interaction [? ]. 

Using Bloch theorem methods to take advantage of the lattice periodicity, the associated eigenvalues 
and eigenvectors of Equation (10) are found to be 

e±(k) = ±t|/(k)| (11) 

and 

where e'"'^'^^^ = and /(k) = e"*^ + 2cos(^) e*^. The bands given by Equation (11) are 
shown plotted along the high symmetry points of the graphene Brillouin zone by the solid lines in 
Figure 3 a. 

We note that since the orbitals in carbon contain one electron each, then the band is half full and the 
Fermi energy of the undoped system is Ep = 0. This is exactly the point where the two bands given by 
e_ and e+ touch and the DOS vanishes. Graphene is thus a zero-bandgap semiconductor, or a semi-metal. 
The resultant Fermi surface consists of six discrete points, only two of which are unique, lying at the 
vertices, or K points, of the hexagonal first Brillouin Zone (BZ) illustrated in Figure 3. Many of the 
interesting properties belonging to graphene arise due to the shape of the bands near the Fermi energy. 
Unlike in conventional semiconductors where the bands are parabolic, the band structure of graphene 
is linear near Ep. This results in electrons or holes near the Dirac points having zero effective mass 
and behaving like relativistic particles which can be described using the Dirac equation from Quantum 
Electrodynamics. It is worth investigating briefly how graphene electrons are approximated in the linear 
regime and the range of energy values over which the approximation is valid. 

The region of the graphene spectrum around Ep = can be approximated linearly to simplify the 
calculation of physical properties. For a point in reciprocal space near the Dirac point, k = K + (5k, 
we find 

e± (|5k|) = ±^!^ |5k| = ±hvp \5k\ (13) 
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Figure 3. (a) shows the band structure of graphene calculated using both the tight-binding 
Hamiltonian (solid lines) and the linearised Hamiltonian (dashed lines). The bands are 
plotted along the high symmetry points of the graphene Brillouin Zone shown by the dashed 
red plot in (b). The shaded rectangular area in (b) is a possible altemative to the standard 
hexagonal Brillouin zone. 

a) b) 





where the Fermi velocity of graphene, vp = \^ ~ 10^ ms"^. In Figure 3a we compare the linear 
approximation of the band structure (dashed lines) with that calculated previously (solid lines). It is seen 
to be a very good approximation in the region surrounding = 0, but loses accuracy quickly outside 
this regime. 

To calculate the coupling using the expressions in the previous section, we need the Green function 
connecting two sites on the pristine graphene lattice. There are a number of different ways to calculate 
this quantity and use it to calculate the coupling, but for illustration we will focus on a real space, energy 
dependent Green function calculated from the dispersion relation and eigenvectors of the non-linearised 
tight-binding Hamiltonian. The Green function we require is given by an expression of the form 



9ir' = (rj.^i l^(^) |r«,n;) 



a a 



2tt iir 



dL, / dk 



(14) 



where 



NiE) 



E if 
t/(k) 

tr(k) 



n 



Hi 



n,: 



1 -^i 
= '2,ni 



2 
1 



and the reciprocal space integral is over the first BZ, where for convenience, we can define a rectangular 
BZ consisting of segments drawn from multiple neighbouring BZs whose area equals that of the 
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hexagonal BZ. A possible rectangular BZ is shown by the shaded area in Figure 3b. We note that 
the Green function expression takes a different form for sites on the same or on opposite sublattices. 
The most obvious, if not most efficient, way to proceed from here is a brute-force two-dimensional 
numerical integration. However, we have shown previously [67] that it is possible to perform one of the 
integrals completely analytically, greatly reducing the numerical cost of the Green function calculation. 
Furthermore, for separations greater than a few lattice spacings in the high symmetry zigzag and armchair 
directions it is possible to approximate the remaining integration to a high degree of accuracy using the 
Stationary Phase Approximation (SPA), to yield a closed form expression for the Green function 



where A{E) is an energy-dependent coefficient and Q{E) can be identified with the Fermi wavevector 
in the direction of separation. The full analytical expressions for these quantities for armchair and zigzag 
separations can be found in [67]. Once calculated using any of these methods, the Green functions can 
be substituted directly into the RKKY integration in Equation (6). Alternatively, they can be used to 
calculate the spin-dependent Green functions required for the non-perturbative total energy expression 
in Equation (4), where the Dyson equation is used to add a spin-dependent potential at the magnetic sites. 

At this stage, it is worth mentioning some of the alternative methods used to calculate the graphene 
Green function for RKKY calculations in the literature. Many of the earlier, analytic studies employ 
the linear approximation to the graphene band structure [57,58,68-70] and may also use time-dependent 
Green functions in lieu of the energy dependent form shown here [57-59]. There is some discrepancy 
between early results for both the sign and decay rate of the interaction (e.g., [57,58,69]), but many of the 
key features to be discussed in the next section are present in [57] and later confirmed by other studies. 
Numerical investigations using Green function techniques have tended to use an energy dependent form, 
but calculated in a variety of ways [67,72? ]. Sherafati and Satpathy [71] employ both full numerical 
integration of the Green function integral in Equation (14) and the recursive Horiguchi method, which 
relates the graphene Green function to those of the triangular lattice, which in turn are expressed in 
terms of elliptic integrals [73]. Although more efficient than the full integration, the Horiguchi method 
was only stable for small values of separation [71]. In the same work, an analytic treatment in the 
linear dispersion regime using energy dependent GFs confirmed, with some refinement, many of the 
features predicted previously using time dependent GFs. In a study of RKKY interactions in disordered 
graphene sheets, Lee et al. [72] used the Kernel Polynomial Method [74] to express the Green functions 
in terms of Chebyshev polynomials with recursively calculated coefficients. Working outside the Green 
function formalism the interaction can be calculated by determining the total energies of the FM and 
AFM alignments of the moments. Such an approach can use tight-binding models similar to those 
presented here [65,66] or more complete ab initio calculations [60-64]. Such methods often rely 
on periodic boundary conditions and relatively small unit cells, leading to problems like interference 
between moments in neighbouring cells [75] and difficulty in calculating the interaction for separations 
greater than a few multiples of the lattice parameter. 



Qd{E) 



(15) 



3. Interaction between Simple Magnetic Impurities 
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In this section we consider localised magnetic moments that are associated with a single lattice 
site in graphene. These can be thought of as substitutional magnetic impurities replacing a single 
carbon atom in the lattice, or as moments induced at a lattice site by the presence of an adsorbate. 
This type of impurity is the most commonly considered in the literature and the exchange interactions 
between pairs of these impurities have been investigated using many of the methods discussed in the 
previous sections. The expressions for the couplings given in Equations (4) and (6) are easily applied to 
moments of this kind as the Green functions required are either those for pristine graphene with a simple 
spin-dependent perturbation at the magnetic sites, or without this perturbation if the RKKY expression 
is used. Throughout this section we will consider two magnetic impurities occupying sites A and B in 
the graphene lattice, separated by a distance D, with an indirect exchange interaction Jab- 

In the earliest studies, the RKKY interaction in graphene was often considered in conjunction 
with Friedel Oscillations (FO). The latter involve distance dependent charge fluctuations relative to an 
impurity site and there is a large degree of overlap in the physics involved. Using an analytic Dirac 
linearisation approach to study FO, Cheianov and Fal'ko predicted that charge density FO would decay 
as 5p ~ D~^, but the RKKY interaction as J ~ [69]. Wunsch et ah, meanwhile, found the 
magnitude of induced moments to decay as D^^ away from a magnetic impurity and predicted the 
RKKY interaction to behave likewise [76]. Before these reports, an early study by Vozmediano et al. 
investigating possible ferromagnetism in graphene layers arising from localised moments near lattice 
distortions had examined the RKKY interaction between moments surrounding elongated cracks in 
graphene [68]. Although more complicated structures than considered elsewhere, a decay rate of 
D"^ for a non-oscillatory, always ferromagnetic interaction was reported. An important paper by 
Saremi clarified many of the features of the interaction [57]. Firstly, a general proof was given 
regarding the sign of the RKKY interaction in a half-filled bipartite lattice. This showed that, due 
to electron-hole symmetry in such systems, the coupling was ferromagnetic between moments on the 
same sublattice and antiferromagnetic between moments on opposite sublattices. In the same work, 
analytic calculations were performed in the linear dispersion regime to determine the functional form 
of the RKKY interaction. These calculations included the explicit use of a cut-off function to prevent 
divergences arising from high energy contributions. There has since been some debate about the form, 
or even the necessity, of the cut-off function and of the nature in which the time-dependent Green 
function integrations are performed [58,59,65,67,71,77? ]. However, many of the principal features 
of the interaction reported by Saremi have been verified by other studies, including a contemporaneous 
paper by Brey et al. [70]. 

1 . FM (AFM) interaction between moments on the same (opposite) sublattice(s), with the magnitude 
of the AFM interaction three times greater; 

2. A decay rate of Jab ~ D~^; 

3. An oscillatory term of the form 1 + cos{2Q{Ef)D), where Q{Ef) is the component of the Fermi 
wavevector parallel to the separation direction. 

Bunder and Lin suggested that differences between the real time and imaginary time approaches 
could lead to a breaking of the theorem regarding the sign of the interaction [58]. In particular they 
report sign changing oscillations for the interaction in zigzag-edged graphene nanoribbons. However, 
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fully-numerical calculations of the coupling by Black-Schaffer do not find such sign changes, which 
the author attributes to the perturbative nature of the previous calculations [65]. In the same work, the 
main predictions made by Saremi are verified using a numerical exact diagonalisation method to find the 
energy difference between the FM and AFM orientations. A minor refinement is the introduction of a tt 
phase shift in the cosine term for moments on opposite sublattices. In fact, Sherafati and Satpathy showed 
that this phase shift is not a constant but depends on the direction of separation, taking a maximum value 
of TT for zigzag separations and vanishing for armchair separations [71]. This additional phase shift term 
can be added to the three features proposed by Saremi to provide a complete description of the RKKY 
interaction between simple single-site magnetic moments in graphene. 

Using the stationary phase form of the Green function in Equation (15), we have previously 
demonstrated how the principal features of the interaction in the RKKY approximation can be 
calculated [67]. The integration over energy in Equation (4) can be reduced to a sum over Matsubara 
frequencies. The functions B{E,e) = A'^{E,e) and Q{E,e) are expanded around Ep and in the low 
temperature limit T — 0, we find 

Jba ~ Im 5^ exp (2zQ(E^) D) (16) 
i 

where 

MEf) - (2Q'(E^))^+i ^^^^ 

is the distance-independent coefficient for the i-th term in the power series, ^ is a non-negative integer 
and B^^\Ef) is the ^th order energy derivative of B{E) evaluated at Ep, resulting from its Taylor 
expansion. In general the leading term in the series should determine the asymptotic decay rate of the 
coupling. For the undoped case it can be shown that the coefficient B^^\Q) = 0, so that the / = 1 term 
dominates and J{Ef = 0) ~ D^'^. The exact oscillatory term predicted by Saremi is reached by forcing 
the real part of the Green function to vanish at Ep = 0. 

Figure 4 shows a numerical calculation of the coupling, calculated using Equation (5), for 
substitutional magnetic impurities in graphene. The coupling is plotted as a function of separation 
for both armchair and zigzag directions, and for impurities on the same and on opposite sublattices. 
Results are shown for both undoped (Ep = 0.0) and doped (Ep = 0.1 |t|) cases. The full range of 
features discussed in this section are present. In particular we note the D^^ decay rate and absence of 
sign-changing oscillations in the undoped cases. We note the presence of a period-3 oscillation in the 
zigzag direction compared to the monotonic decay for the armchair case. This results from slightly 
different commensurability effects in the two directions. In undoped graphene, the oscillations are 
determined by cos^(kF ■ D). The components of the Fermi wavevector in the armchair and zigzag 
directions are and |^ respectively. Since the spacing between unit cells is ^/Sa for armchair and 



a 



for zigzag separations, the cosine squared term takes a constant value of 1 in the armchair direction but 



cycles repeatedly through the sequence |, |, 1 as the separation is increased in the zigzag direction. Also 
evident in the zigzag separation plot is the vr phase shift between same and opposite sublattice cases. 
When the system is doped, as in the bottom panels, a longer-ranged interaction decaying as J ~ D^^ is 
seen. Long-wavelength, sign-changing oscillations are also introduced by the change in the shape of the 
doped Fermi surface. 
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Figure 4. The indirect exchange coupling between two moments separated in the armchair 
(direction) are shown in the left (right) panel as a function of moment separation. Results for 
the undoped (doped) case are shown in the top (bottom) panels. Results are shown for both 
the same sublattice (black, solid lines) and opposite sublattice (red, dashed lines) cases. 
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4. Beyond the Simple Case 

In the previous section we highlighted the features of the magnetic interaction in graphene for 
the simplest kind of impurities. These impurities are the kind usually investigated within the RKKY 
approximation. In this section we consider a number of extensions to the simple Hamiltonian describing 
the system and examine the effects they have on the coupling. We consider first a more general or 
realistic description of the magnetic impurities and examine the effect of impurity parameterisation on 
the interaction. We show that the use of the pristine Green function in the RKKY approximation ignores 
these effects which can in principle overturn some of the results noted in the previous section. Then 
we move to the case when the impurity connects to more than one atom in the graphene lattice and 
examine the effect on the decay rate of the interaction. Finally, more detailed descriptions of the graphene 
electronic structure are considered. We discuss the possible effects of electron-electron interactions and 
spin- orbit coupling and examine the results of some ab initio calculations. 



4.1. Impurity Parameterisation 

Calculations using the RKKY approximation do not generally account properly for the local spin 
dependent potentials describing the magnetic moments. In fact, the parameterisation of the magnetic 
impurity is generally neglected through the usage of the pristine Green functions. In this section 
we consider the effect of three parameters that can be used to characterise the magnetic impurities. 
In addition to the magnetic moment (m) and band-centre shift (5), shown in Figure 2 and which 
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can be calculated, for example, using self-consistent mean-field calculations, we also consider the 
hopping potential between the lattice carbon sites and the impurity site {t') which should differ from the 
carbon-carbon hopping. This approach is similar to the Anderson model describing localised magnetic 
impurity states in metals [78]. These parameters will vary between the different magnetic species that 
can be chosen as the embedded impurities. 

In Figure 5 we plot the coupling between substitutional impurities on the same sublattice, calculated 
using the Equation (4) and exact Green functions, as a function of separation along the armchair 
direction for three different parameter sets (m, 5, t'). The first of these (0.6, 0.0, t) closely replicates 
the results of the RKKY approach as it considers only a band-splitting, has no band-centre shift and 
uses the carbon-carbon hopping value. The numerical calculations in the previous section used a similar 
parameterisation. The middle and bottom plots use the parameters (0.6, 5.0, 0.8t) and (0.6, 8.0, 0.6t) 
respectively. In these plots we see the formation of an unusual feature not predicted by the RKKY 
approximation. For quite a large range of distances we note a preferential antiferromagnetic alignment 
between the moments before the sign flips and the standard ferromagnetic coupling with a decay is 
recovered. A similar sign-changing behaviour has been reported in nanotubes [49] and also by ab initio 
calculations attempting to probe the interaction in graphene [64]. It is worth examining further how this 
feature depends on the parameterisation of the magnetic moments. 

Figure 5. The magnetic coupling between two magnetic sites on the same sublattice as 
a function of their separation, D, along the armchair direction for three different impurity 
parameterisations: (a) (m = 0.6, 5 = 0.0, t' = t); (b) (m = 0.6, 5 = 5.0t, t' = O.St) and 
(c) (m = 0.6, 6 = 8.0t,t' = O.Qt). The insets show log-log plots where a sign change in the 
coupling is evident from a dip feature. 
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In Figure 6 we present a number of phase diagrams showing the sign and strength of the coupling 
for different values of these parameters. Each diagram represents an area of (m, S) phase space with 
ferromagnetic (antiferromagnetic) couplings given by a blue (red) shading that is darker for larger 
magnitude couplings. The diagrams on the top row correspond to a separation of 10^/a between the 
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magnetic moments and from left to right show the cases of t' = l.Ot, 0.8t, 0.4t. The bottom panels 
show the same cases for a larger separation of AQy/a. By examining the border between the blue and 
red regions in these plots we can infer under what circumstances the sign-change behaviour described 
above occurs. In all cases the border position varies only weakly with the magnetic moment (m) 
or hence the band- splitting {Vex)- A stronger dependence is found on the band-centre shift {5) and 
we find that, in general, an anti-ferromagnetic alignment is found above a critical value of 5. The 
band-centre shift is strongly dependent on the occupation of the magnetic orbital. For a bipartite 
lattice like graphene a substitutional impurity with a half-filled orbital gives zero band-centre shift 
when t' = t. As we move away from half-filling a larger band-centre shift is required to return the 
correct band occupation. For smaller values of t' we note that the border between the blue and red 
regions shifts towards the left, meaning that smaller band-centre shifts will lead to an AFM alignment. 
Increasing the distance between the impurities reduces the phase-space area corresponding to an AFM 
alignment by shifting the border to the right and requiring larger band-centre shifts. This finding agrees 
with the distance-dependent plots in Figure 5 for fixed parameters which show that in the asymptotic 
limit the coupling changes sign and returns the FM alignment predicted by the RKKY approximation. 
However, as in the bottom panel of Figure 5, the magnitude of the coupling has essentially decayed 
to zero before the sign change occurs so the only significant coupling between two such impurities is 
antiferromagnetic. At this point, we note that the theorem proposed by Saremi [57] suggesting that the 
same (opposite) sublattice coupling is always FM (AFM) only applies for the case of bipartite lattices 
with half-filling. The configurations discussed here break the half-filling requirement and electron-hole 
symmetry, at least locally, and this has a dramatic effect on the nature of the interaction. In particular, 
depending on the moment parameterisation, a strong antiparallel alignment between two substitutional 
magnetic impurities in graphene may persist to considerable separations. This result is contrary to many 
predictions made regarding the RKKY interaction in graphene and shows that caution needs to be applied 
when extrapolating the results of simple calculations to more realistic magnetic impurities. 

4.2. Multi-Site Impurities 

The features of the RKKY interaction discussed in Section 3 contain a strong sublattice dependence. 
A substitutional or top-adsorbed impurity is associated with a single site on the graphene lattice and 
therefore also with a single sublattice, as illustrated in Figure 7a. We have seen how the sign of the 
interaction, at least for simple impurity parameterisations, depends on whether the interaction is between 
sites on the same or on opposite sublattices. However, another type of impurity can be considered which 
is associated with an equal number of sites on both sublattices. Two examples of this type of impurity are 
shown in Figure 7. These are the bridge impurity (b), adsorbed onto the lattice midway between an atom 
from each sublattice, and the plaquette, or centre-adsorbed, impurity at the centre of hexagon (c) and 
bonding equally to the six surrounding carbon atoms. Plaquette impurities have been investigated more 
thoroughly than bridge impurities and a number of interesting predictions have been made. An important 
difference is noted between the results for incoherent and coherent type moments [57]. For incoherent 
moments, the calculation essentially consists of adding the contributions of the thirty-six single-site 
interactions between the associated sites of one moment and the other. Coherent moment calculations 
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Figure 6. (m, 5) phase-space diagrams for separations of 10-\/a (top row) and 40-\/a (bottom 
row) for three different values of the impurity-carbon hopping parameter t' . The sign of the 
coupling is indicated by the colour (blue for FM and red for AFM) and the strength of the 
coupling by the degree of shading. 
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are more complete and can be performed, for example, using Equation (4) with the spin-dependent Green 
functions of a graphene sheet with the two impurities connected to the graphene sheet using the Dyson 
equation and an appropriate perturbation potential. Saremi noted a decay rate of for incoherent 
plaquette impurities [57], a result confirmed by further studies [65,70,71]. However, Saremi also noted 
that for the coherent case the D^^ term vanishes and that the decay in this instance is faster. This 
is similar to the decrease in the interaction range reported for plaquette impurities in metallic carbon 
nanotubes [48], where the decay rate was found to be compared to the rate predicted for 
substitutional or top-adsorbed impurities. Uchoa et al. have since reported a decay rate of D^^ for 
plaquette impurities in graphene [79]. 

Figure 7. The three types of impurity configuration discussed: (a) substitutional atom; 
(b) bridge- adsorbed atom and (c) plaquette impurity. 

(a) (b) (c) 
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The sign of the interaction between plaquette impurities is also of interest. An always AFM coupling 
is widely reported in the literature [57,65,71], although Uchoa et al. report a FM coupling for smaller 
values of separation [79]. The prevalence of an AFM coupling can be understood from the single- site 
coupling results. We may initially expect the plaquette impurity to average out the sublattice-driven 
preference for an always FM or AFM coupling. However, we recall that the magnitude of the AFM 
coupling in the single site case is three times larger than that of the FM coupling. Thus the AFM 
contributions from inter- sublattice terms are larger than FM contributions from intra- sublattice terms 
and an overall AFM coupling is found between plaquette impurities. Another interesting feature for this 
impurity type is that the period-3 oscillations noted for zigzag separations are no longer present and a 
smooth monotonic decay of the coupling is noted for separations in all directions. The jagged features 
are averaged out by the fact that a plaquette impurity is associated with an equal number of sites from 
each group defined by the period-3 sequence noted for single- site impurities. 

Sherafati and Satpathy [71] have examined the case of bridge impurities and find a rare example of 
sign-changing oscillations in the interaction in undoped graphene. They report that every third value of 
separation has an AFM coupling, and the other two are FM. 

4.3. Spin-Orbit Effects, Electron-Electron Interactions and ab initio Calculations 

Having spent some time examining the effects of more detailed impurity parameterisations, it is now 
worth considering the role a more complicated description of the graphene electronic structure may 
play. Such a description may include additional terms in a model Hamiltonian calculation or a full 
ab initio treatment using a Density Functional Theory (DFT) approach. An early study of this type was 
undertaken by Dugaev et al. [80] and investigated the interaction in a graphene system where a small 
gap was introduced by the inclusion of spin-orbit effects. Using an analytic approach based on the Dirac 
equation, the prediction of an always FM interaction with a long tail exponential decay was made. 

Black-Schaffer [66] included the effect of electron-electron interactions in a numerical study of 
the RKKY interaction. These interactions were included using a mean-field solution of a Hubbard 
Hamiltonian with a finite onsite Coulomb repulsion term U at each site in the graphene lattice. Similar 
methods have been used previously to investigate intrinsic magnetism in graphene arising from the spin 
polarisation of localised states near zigzag edges. A number of interesting effects are noticed as the value 
of U is increased. The oscillatory features in the zigzag direction are removed, as is the factor of three 
amplitude difference between inter- and intra- sublattice couplings. Furthermore, the decay exponent 
decreases from 3 for j = to below 2 when j > 2. At this value of U the interaction is also isotropic, 
with no difference between the interactions for armchair and zigzag separations. 

A few ab initio studies have been performed which address exchange couplings between localised 
magnetic moments in graphene. These studies are generally difficult to compare with the other results 
presented here for a number of reasons. Due to a much larger computational cost, ab initio calculations 
are generally only performed for very small values of separation. The majority of these calculations 
use periodic boundary conditions, which can result in magnetic interactions between impurities in 
neighbouring unit cells [75,81]. In fact, such interactions can lead to unusual results in ab initio studies 
examining a single magnetic impurity due to the non-oscillatory nature of the RKKY interaction in 
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graphene. Preferential AFM alignment between impurities in neighbouring unit cells can potentially 
lead to the total suppression of the magnetic moment [31,75]. Periodic unit cell calculations within a 
tight-binding scheme can minimise this effect by including a large "padding" region around the 
impurities [65], but this is usually not feasible in DFT. Signatures of the RKKY interaction can also be 
seen in DFT calculations when the effective distance between moments in single impurity concentrations 
is varied by increasing the size of the unit cell [32,60]. From a study of pairs of localised moments 
arising from simple vacancy defects, Pisani et al. [61] extracted a decay rate of D~^ '^^ for moment 
separations of up to 25A. This agrees with Black-Schaffer's calculations including electron-electron 
interactions when the onsite Coulomb repulsion term is j = 2.1. A study by Santos et al. [62] 
however extracts a decay rate of D^"^ "^^ for a FM decay between substitutional cobalt impurities on the 
same sublattice. This result is significantly closer to the rate predicted without electron- electron 
interactions. Calculations between Fe atoms meanwhile suggest a preferred AFM alignment between 
impurities on the same sublattice [75], suggesting that the nature of the interaction depends strongly 
on the magnetic species chosen, in agreement with the findings of Section 4 of the current work. 
Moments arising in hydrogenated armchair-edged graphene nanoribbons generally preferred an AFM 
alignment [63]. Carbon adatoms, adsorbed in the bridge configuration, were found for small separations 
to have a FM alignment. For larger separations both FM and AFM alignments were found for different 
values of separation [64]. 

These findings show that there is still quite a large amount of discrepancy between the interactions 
predicted using simple RKKY-like models and more complete DFT calculations, which appear very 
strongly dependent on impurity configuration. Furthermore, direct comparison of results from both 
methods is often thwarted by the separation constraints imposed by DFT calculations and the possibly 
oversimplistic treatment of magnetic impurities in model calculations. 

5. Manipulating the Interaction 

In this section we discuss a number of factors which may allow manipulation of magnetic interactions 
in graphene. A degree of control over the properties of the interaction may be required for spintronic 
applications, for example to switch on and off spin currents or to dynamically change the magnetic 
ordering in a system. Another important consideration is that the interactions discussed so far have been 
very short-ranged, due to a fast decay rate arising from the peculiar electronic dispersion relation in 
graphene. Methods of augmenting the interaction, or slowing the rate of decay, may thus help overcome 
some of the difficulties in realising magnetically-doped graphene devices. 

5.1. Geometry Effects 

A number of studies have looked at magnetic interactions in graphene systems with geometries 
other than the standard two-dimensional sheets considered so far in this work. These include 
nanotubes [47-50,82,83], nanoribbons [58,63,65,66,84,85], nanoflakes [87] and multilayer graphene 
systems [59,88-91]. Quasi-one-dimensional graphene systems, such as nanotubes and nanoribbons, 
are expected to display longer ranged interactions than those present in two-dimensional graphene 
sheets. Indeed, metallic nanotubes have been shown to display a long ranged D^^ interaction when 
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substitutional, top-adsorbed or bridge-adsorbed impurities are considered [47,48]. However, similar to 
the case of graphene, the interaction between plaquette impurities is found to decay at the much faster 
rate of D~^. The sublattice dependent sign rules discussed for graphene also hold for nanotubes within 
the RKKY approximation, but can be broken when the parameterisation of the impurities is altered [49]. 
The eigenstates of carbon nanotubes are subject to periodic boundary conditions in the circumferential 
direction which ensure the equivalence of all lattice sites. This is not the case in nanoribbon systems 
where edges are present and each lattice site can be characterised by its distance from the edges. Thus 
we should expect the interaction to contain a dependence on the location of the impurities across the 
width of the ribbon. Such a dependence is found in many properties of doped ribbons, including binding 
energies [92], conductance [93,94] and the magnitude of magnetic moments [95]. The presence of 
localised edge states for zigzag edge geometries [17-19] should also be expected to introduce features not 
present in graphene sheets or nanotubes. Black-Schaffer [65] examined the interaction in nanoribbons, 
finding that armchair edges did not significantly alter the coupling from the bulk case. However, for 
localised moments at the edge of zigzag ribbons, an exponentially decaying interaction was reported. 
The interaction was FM for moments along one edge and AFM for opposite edges. However, the 
magnitudes of the FM and AFM interactions were equal, unlike the larger AFM magnitude found for the 
bulk case. Towards the ribbon centre, the results for the bulk were recovered. Including electron-electron 
interactions in the calculation suggested that the distance dependence vanished for reasonable values of 
U, leading to an extremely long-ranged interaction [66]. Szalowski studied the RKKY interaction in 
finite graphene nanofiakes, finding that for certain geometries the addition of a single additional charge 
carrier could change the sign of the interaction [87]. A recent work by Klinovaja and Loss extends the 
study of the interaction in one-dimensional graphenes to include spin-orbit interactions[86]. The RKKY 
interaction has been studied in bilayer graphene systems by several authors [59,88-90]. The quadratic 
nature of the dispersion relation in bilayers means that many of the features present in monolayer 
graphene are removed. Another interesting feature is that, for A-B stacking, a plaquette impurity in 
one layer is embedded above a single graphene lattice site on the other layer. Extending the study to 
more than two layers, Jiang et al. [91] find a number of different decay rate possibilities, depending on 
the positions of the two impurities and on the number of layers. 

5.2. Doping and Disorder 

The key features of the RKKY interaction when the Fermi energy is shifted away from half-filling 
are clear from the bottom panels of Figure 4. These are namely a longer ranged interaction that decays 
as D^"^ and sign-changing oscillations, the periods of which depend on the component of the Fermi 
wavevector in the direction of separation. The appearance of oscillations is due to the breaking of the 
commensurability effect between the Fermi wavevector and the lattice spacings. Many of the studies 
of the lEC in graphene establish these key features in doped or gated graphene, but they have been 
explored in more detail by Sherafati and Satpathy [96], where the possibility of controlling the sign of 
the interaction with a gate voltage is raised. The possibility of controlling the interaction using a gate 
voltage has also been discussed for bilayer graphene [89]. 



Crystals 2013, 3 



21 



Lee et al. [72] examine the effects of nonmagnetic disorder on the RKKY interaction in monolayer 
graphene. The Anderson model was employed to study both onsite and hopping term disorder. Onsite 
disorder was found to induce sign-changing oscillations in the coupling in certain directions, whereas 
off-diagonal disorder only effected the amplitude, and not the sign, of the interaction. This suggests that 
sign-changing behaviour is introduced by the breaking of sublattice symmetry, in a similar manner to that 
discussed for impurity parameterisation in Section 4. For weak disorder, the interaction retained similar 
decay behaviour to the clean case. In the strongly disordered, localised regime exponential suppression 
of the coupling with increasing moment separation was observed. In a recent work, the same authors 
examine the interplay of disorder and gate voltage [97]. 

5.3. Strained Graphene Systems 

Strain-engineering of graphene systems has received a lot of attention in recent literature due to the 
possibility of tuning many of the physical properties of graphene [98-109]. The degree of tuning is 
enhanced by the different types of strain that can be applied. Apart from simple uniaxial strains [98,1 10], 
more exotic features like creases and bubbles can be introduced [101,104,106,111-113]. In a recent 
work, we have explored the possibility of manipulating the lEC in graphene by the application of uniaxial 
strains parallel or perpendicular to the moment separation direction [109]. Another recent work also 
explores this topic [108]. We find that both amplification and suppression of the magnetic coupling 
can be achieved. We reported a general trend of amplification for strain perpendicular to the moment 
separation direction, or suppression for strain parallel to this direction. Also noted are oscillations in the 
amplification as strain is increased for moments separated in non-armchair directions. Such oscillations 
suggest the intriguing possibility of selectively turning on or off the coupling between moments and also 
of controlling the inter- and intra- sublattice couplings independently. Since the lEC underpins physical 
features, including overall moment formation and magnetotransport response, the ability to fine tune the 
coupling with strain may lead to interesting spintronic applications. Thorough investigation of graphene 
systems with magnetic impurities and with different types of strain applied may yield a further range of 
tuneable spintronic properties. 

5.4. Spin Precession and Dynamic RKKY 

One of the major obstacles in the path of exploiting the RKKY interaction between magnetic 
impurities in graphene for spintronic devices is the very short range of the interaction. It has been 
suggested that the range of RKKY-like interactions can be augmented if the magnetic moments are 
set in motion, for example, if they are set to precess by the application of a suitable time-dependent 
magnetic field [114,115]. Such an interaction is driven by non-equilibrium spin currents emanating 
from the precessing moments [116]. The magnitude of the interaction can be measured by a quantity 
called the dynamic spin susceptibility, which describes the response of the magnetism of the system to 
a dynamic magnetic perturbation. Investigation of this quantity in carbon nanotubes [82] has revealed a 
decay rate slower than that of the static coupling. Further studies of spin dynamics in graphene systems 
have suggested the use of these materials as spin waveguides [117], spin-pumping transistors [118] and 
spin current lenses [119]. Recent experimental evidence also suggests possible long-range spin current 
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behaviour in graphene [120]. In a recent study [121], we examined the spin susceptibility in graphene as 
a dynamic analogue of the static RKKY coupling, with a particular focus on the separation dependence 
of the interaction. We find that the decay rate for the dynamic RKKY approaches for both 
doped and undoped systems, giving a significant augmentation of the interaction range. Furthermore, 
the dynamic interaction can be related to fluctuations in the lifetimes of magnetic excitations at the 
impurity sites. Such quantities have been probed in other materials using inelastic scanning tunnelling 
spectroscopy [122-125], suggesting that this method may be used to detect signatures of the dynamic 
RKKY interaction in graphene. 

6. Conclusion and Experimental Considerations 

In this review we have examined many aspects of indirect exchange interactions in graphene 
systems. The various approaches to calculating the coupling, ranging from analytical solutions of 
the RKKY approximation using the Dirac approximation, to numerically expensive ab initio total 
energy calculations, via tight-binding energy difference calculations with varying levels of impurity 
parameterisation, were discussed. Each of these methodologies has merits and pitfalls, and the 
contrasting results suggest that there is still some work to be done before useful comparison can be made 
with experimental results. On the experimental side, it is understandable that with a decay rate as fast 
as it is difficult to probe the interaction for any reasonable separation. The presence of magnetism 
in disordered graphene systems may indicate the presence of an exchange coupling between magnetic 
moments formed around defects. Nuclear magnetic resonance experiments reveal that these defects 
have indeed magnetic moments, since they couple to implanted Fe atoms [126]. However, whether or 
not these moments couple with each other, or with the graphene lattice, to form a ferromagnetic state 
is a controversial subject and many of the results in this area have proved difficult to reproduce [5]. 
Progress in spin-dependent scanning tunnelling spectroscopy has enabled the RKKY interaction to 
be measured directly from the magnetisation curves of individual magnetic adatoms [127,128]. As 
discussed above, signatures of the dynamic analogue of the RKKY interaction may be detectable from 
the magnetic excitation lifetimes of individual moments measured using inelastic scanning tunnelling 
spectroscopy [122-125]. In this case, the longer range of the interaction may make experimental 
detection more feasible [121]. Using a gate voltage to increase the interaction range may also be 
necessary to probe the interaction in the static case. To our knowledge these types of experiment have 
not yet been performed with a graphene host medium. Signatures of indirect exchange interactions 
may also be present in transport measurements taken in magnetically-doped graphene systems [49,129], 
since different orientations of the magnetic moments lead to different scattering regimes for up- and 
down- spin electrons. 

Future work on the topic of indirect exchange interactions in magnetically doped graphene should 
focus on the gaps and discrepancies between the various theoretical models, and between these 
models and possible experimental configurations. In particular, the strength of the interaction for 
specific magnetic impurities and whether it survives from the ideal model to experimentally realisable 
conditions of finite temperature and lattice disorder are important considerations. Another area of 
exploration is the interplay between magnetic impurities and other features that can be introduced into 
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the graphene lattice. The topics of different geometries, in the form of nanoribbons and nano tubes, and 
of geometrical distortions, in the form of strain, were briefly discussed in this review and results were 
mentioned for simple cases. More complicated geometries or distortions may introduce new features 
and more possibilities to manipulate the interaction. In moving towards technological application of 
RKKY-like interactions, the study of larger ensembles of magnetic impurities will be necessary. The 
magnetic and electronic properties that emerge from such ensembles should be dictated by the pairwise 
indirect exchange interactions discussed in this review, and may lead to a range of interesting features 
and applications. 
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